classdef edffile < handle
   
    methods
        
        function self = edffile(path)
            fh = fopen(path,'r','n','US-ASCII');
            self.fs = fread(fh,inf,'char=>char')';
            fclose(fh);
        end
        
        function co = coordinates(self)
            % Parse the file
            % see
            % http://stackoverflow.com/questions/747296/how-to-elegantly-ignore-some-return-values-of-a-matlab-function
            if isempty(self.coo)
                [ans,ans,ans,ans,coa] = regexp(self.fs,['([0-9\.])+\s+([0-9\.]+)'...
                    '\s+([0-9\.]+)\s+([0-9\.]+)\s+\.{3}'],'dotexceptnewline');
                co = zeros(numel(coa),4);
                for i = 1:numel(coa)
                    co(i,:) = str2double(coa{i});
                end
                self.coo = co;
            else
                co = self.coo;
            end
        end
        
        function s = signal(self)
            % Gets the signal with out saccades
            
            s = coordinates(self);
            saccades = self.getSaccades();
            p0 = s(1,1)-1;
            saccades = saccades - p0;
            s = s(:,2);
            % Reports saccades as non values
            for i = 1:size(saccades,1)
                s(saccades(i,1):saccades(i,2)) = NaN;
            end
            p1 = start(self);
            p2 = final(self);
            s = s(p1{1}{1}-p0:p2{1}{1}-p0);
            
        end
        
        function v = velocity(self,signal,cutoff,sampling)
            % Filter the data at 20 Hz and takes the gradient
            if nargin < 3; cutoff = 20; end
            if nargin < 4; sampling = 1000; end
            
            [v,b,a] = self.filter(signal,cutoff,sampling);
            self.transferFunction.zeros = b;
            self.transferFunction.poles = a;
            v = gradient(v(1:end));
            v = v(:); 
        end
        
        function saccades = getSaccades(self)
           % Returns a n by 2 matrix, where n is the total number of saccades 
           % in the first column is the initial time and in the second
           % column is the end time.
           [ans,ans,ans,ans,sacc] = regexp(self.fs,...
               'ESACC\s+[RL]\s+([0-9\.]+)\s+([0-9\.]+)','dotexceptnewline');
           saccades = ones(numel(sacc),2);
           for i = 1: numel(sacc)
              saccades(i,1) = str2double(sacc{i}{1});
              saccades(i,2) = str2double(sacc{i}{2});
           end
        end
        
        
        function msg = messages(self)
            % Finds all the messages
            [ans,ans,ans,ans,msg] = regexp(self.fs,'MSG\s([0-9\.]+)\s+(.+)','dotexceptnewline');
            for i = 1:numel(msg)
                msg{i}{1} = str2double(msg{i}{1});
            end
        end
        
        function h = header(self)
            [ans,ans,ans,h] = regexp(self.fs,'\*\* .+','dotexceptnewline');
            
        end
        
        function s = start(self)
            [ans,ans,ans,ans,s] = regexp(self.fs,'MSG\s([0-9\.]+)\s+START','dotexceptnewline');
            for i = 1:numel(s)
                s{i}{1} = str2double(s{i}{1});
            end            
        end
        
        function e = final(self)
            [ans,ans,ans,ans,e] = regexp(self.fs,'MSG\s([0-9\.]+)\s+END','dotexceptnewline');
            for i = 1:numel(e)
                e{i}{1} = str2double(e{i}{1});
            end            
        end
        
        function hp = plotEyes(self)
            f = self.final();
            s = self.start();
            c = self.coordinates();
            hp = plot(c(:,1)-s{1}{1},c(:,2),'k','linewidth',1); hold on
            title('x-coordinates')
            xlabel('Seconds')
            text(0,c(s{1}{1}-c(1,1)+1,2),'Start');
            scatter(0,c(s{1}{1}-c(1,1)+1,2),'r'); hold on
            if f{1}{1}-c(1,1) + 1 > size(c,1)
                text(f{1}{1}-s{1}{1},c(end,2),'End');
                scatter(f{1}{1}-s{1}{1},c(end,2),'r'); hold on
            else
                text(f{1}{1}-s{1}{1},c(f{1}{1}-c(1,1) + 1,2),'End');                
                scatter(f{1}{1}-s{1}{1},c(f{1}{1}-c(1,1) + 1,2),'r'); hold on
            end
            xlim([-(f{1}{1}-s{1}{1})*0.1,(f{1}{1}-s{1}{1})*1.1])
        end
        
        
        function p = parametersReal(self)
            parameters = {
                'BLOCK'
                'TRIAL'
                'SUBJECT'
                'EXPERIMENT'
                'FPC'
                'OCC_X'
                'OCC_W'
                'JUMP_FRAME'
                'JUMP_dPHASE'
                'FIXATION'
                };
            p = struct;
            for i = 1:numel(parameters)
                [a,b,c,d,e] = regexp(self.fs,[parameters{i} ' (.+)'],'dotexceptnewline');
                p.(parameters{i}) = str2double(e{1});
            end                     
        end
        
    end
    methods(Static)
       
        function signal = interpolateSignal(signal,method)
            if nargin < 2; method = 'spline'; end
            v_nan = isnan(signal);
            x     = 1:numel(signal);
            c = interp1(x(~v_nan),signal(~v_nan),x(v_nan),method);
            signal(v_nan) = c(:);
        end
        
        function [ns b,a] = filter(signal,cutoff,sampling)
            
            % Butterworth filter
            % Implements a Butterworth filter of first order.
            
            % Returns the normalized values of the transferred function
            % See http://www.rowetel.com/blog/?p=2042
            % See http://www.apicsllc.com/apics/Sr_3/Sr_3.htm
            
            % Returs the normalized values of the transferred function of
            % a finete impulse response filter. Uses a blackman window:
            % http://en.wikipedia.org/wiki/Window_function
            n = round((4*sampling)/cutoff);
            a = 1; % FIR filter
            b = fir1(n,(cutoff*2)/sampling,'low',blackman(n+1));
            
            ns = filtfilt(b,a,signal);
            
            ns = ns(:);
        end
        
        function phase = phaseMLE(signal,frequency,sampling)
            % MLE Phase, frequency is the assumed frequency in HZ and 
            % sampling is the sampling rate in herz.
            f = 1/(frequency*sampling);
            % Transform signal into complex signal
            signal_c = hilbert(signal);
            signal_c = signal_c(:);
            % Assuming that the freuncy is known!!!!
            % Sum!
            exp_part = exp(-i*2*pi*f*(0:(numel(signal_c)-1)));
            exp_part = sum(signal_c.*exp_part(:));
            %MLE!
            phase = atan(imag(exp_part)/real(exp_part));
            
        end
         
    end
    
    
    properties 
        fs % File string 
        coo
        % The poles and zeros of the transfer function of the filter
        transferFunction = struct('zeros',[],'poles',[])
    end
    
end